Evolutionary Biology

This paper is concerned with rank deficiency in phenotypic covariance matrices: first to establish it is a problem by measuring it, and then proposing methods to treat for it. Significant rank deficiency can mislead current measures of whole-shape phenotypic integration, because they rely on eigenvalues of the covariance matrix, and highly rank deficient matrices will have a large percentage of meaningless eigenvalues. This paper has three goals. The first is to examine a typical geometric morphometric data set and establish that its covariance matrix is rank deficient. We employ the concept of information, or Shannon, entropy to demonstrate that a sample of dire wolf jaws is highly rank deficient. The different sources of rank deficiency are identified, and include the Generalized Procrustes analysis itself, use of the correlation matrix, insufficient sample size, and phenotypic covariance. Only the last of these is of biological interest. Our second goal is to examine a test case where a change in integration is known, allowing us to document how rank deficiency affects two measures of whole shape integration (eigenvalue standard deviation and standardized generalized variance). This test case utilizes the dire wolf data set from Part 1, and introduces another population that is 5000 years older. Modularity models are generated and tested for both populations, showing that one population is more integrated than the other. We demonstrate that eigenvalue variance characterizes the integration change incorrectly, while the standardized generalized variance lacks sensitivity. Both metrics are impacted by the inclusion of many small eigenvalues arising from rank deficiency of the covariance matrix. We propose a modification of the standardized generalized variance, again based on information entropy, that considers only the eigenvalues carrying non-redundant information. We demonstrate that this metric is successful in identifying the integration change in the test case. The third goal of this paper is to generalize the new metric to the case of arbitrary sample size. This is done by normalizing the new metric to the amount of information present in a permuted covariance matrix. We term the resulting metric the relative dispersion, and it is sample size corrected. As a proof of concept we us the new metric to compare the dire wolf data set from the first part of this paper to a third data set comprising jaws of Smilodon fatalis. We demonstrate that the Smilodon jaw is much more integrated than the dire wolf jaw. Finally, this information entropy-based measures of integration allows comparison of whole shape integration in dense semilandmark environments, allowing characterization of the information content of any given shape, a quantity we term latent dispersion.

Quantitative analyses of morphological variation using geometric morphometrics are often performed on 2D photos of 3D structures. It is generally assumed that the error due to the flattening of the third dimension is negligible. However, despite hundreds of 2D studies, few have actually tested this assumption and none has done it on large animals, such as those typically classified as megafauna. We explore this issue in living equids, focusing on ventral cranial variation at both micro- and macro-evolutionary levels. By comparing 2D and 3D data, we found that size is well approximated, whereas shape is more strongly impacted by 2D inaccuracies, as it is especially evident in intra-specific analyses. The 2D approximation improves when shape differences are larger, as in macroevolution, but even at this level precise inter-individual similarity relationships are altered. Despite this, main patterns of sex, species and allometric variation in 2D were the same as in 3D, thus suggesting that 2D may be a source of noise that does not mask the main signal in the data. However, the problem is complex and any generalization premature. Morphometricians should therefore test the appropriateness of 2D using preliminary investigations in relation to the specific study questions in their own samples. We discuss whether this might be feasible using a reduced landmark configuration and smaller samples, which would save time and money. In an exploratory analysis, we found that in equids results seem robust to sampling, but become less precise and, with fewer landmarks, may slightly overestimate 2D inaccuracies.

For many species, sexual dimorphism is one of the major sources of intraspecific variation. This is the case in some extant great apes, such as gorillas and orangutans, and to a lesser degree in humans, chimpanzees and bonobos. This variation has been well documented in various aspects of these species skeletal anatomy, including differences in the size and shape of the body, cranium, canines, and cresting of males and females, but less is known about sexually dimorphic variation of great ape mandibles. This is particularly important for building robust analog models to interpreting variation in the early hominin fossil record which preserves a large proportion of isolated mandibles and partial mandibles. Here we describe the phenotypical expression of sexual dimorphism in the mandible of six extant hominoid species, including humans, using geometric morphometrics. Our analyses show that the extent of sexual dimorphism in mandibular size and shape amongst the species studied is not the same, as well as the presence of significant differences in the degree of sexual dimorphism being expressed at different sections of the mandible. Furthermore, we find significant differences in how sexual dimorphism is expressed phenotypically even amongst closely related species with small divergence times. We discuss the potential pathways leading to such variation and the implications for extinct hominin variability.

Dental morphology and wear patterns provide insight into the dietary adaptations and ecological niches of living and extinct herbivores. Traditional classification statistics such as Linear Discriminant Analysis (LDA) are limited by assumptions of linearity, normality, and homoscedasticity. This study quantifies mesowear, the shape of molar cusps resulting from occlusal wear, and evaluates the performance of non-linear machine learning models in predicting herbivore diets based on geometric morphometric (GMM) data from adult mandibular second molars (M2) in bovids. We applied Generalized Procrustes Analysis and Principal Component Analysis (PCA) to digitized occlusal shape coordinates from 132 M2 specimens across 64 species. Using the resulting principal component scores, we compared the classification accuracy of LDA with three non-linear models: Random Forest, K-Nearest Neighbors, and Gradient Boosting. While LDA achieved a cross-validated accuracy of just 31%, all non-linear models achieved 99% cross-validation accuracy and 90% test accuracy, demonstrating substantially improved performance. Misclassification analyses revealed that non-linear models more effectively captured complex shape differences, particularly among species with overlapping wear patterns. Our findings support the integration of machine learning with geometric morphometrics to quantify mesowear and improve dietary classification, providing a framework for robust paleoecological inference.

Writers from Darwin to Dawkins have noted that, in humans, women are considered the "beautiful sex," whereas in most species, it is the males who display more elaborate and visually striking traits. This reversal of typical sex roles is a peculiarity of the human species and has been the focus of extensive theoretical debate. Yet, it has never been systematically examined or empirically verified. Here, we present a comprehensive cross-cultural meta-analysis of same-sex and opposite-sex ratings of facial attractiveness from around the globe. Our findings confirm the existence of a robust "Gender Attractiveness Gap" (GAP), with female faces rated significantly more attractive than male faces across rater genders, cultural backgrounds, and portrayed ethnicities. Notably, the effect is more pronounced among female than male raters, suggesting gender-specific modulation. We show that approximately two-thirds of the GAP is mediated by structural facial sex-typicality. However, this mediation is asymmetric: controlling for facial dimorphism substantially reduces the attractiveness of female--but not male--faces, indicating a specific aesthetic preference for structural femininity by both male and female raters. This suggests that attractiveness judgments are driven by general aesthetic evaluation processes that extend beyond heterosexual mate choice. We also observe a general tendency for greater stringency among male raters. Overall, these findings contribute to evolutionary psychology and social perception, offering new insights into mate selection theories and underscoring the importance of accounting for gender-specific and cultural influences in attractiveness judgments.

ObjectivesDental variation within populations and, even more so, within individuals is far less well understood than variation between populations. This is problematic as a single tooth type is often used as a representative of the whole dentition, despite a lack of understanding of intra-tooth type relationships. This research investigates the variation of dental tissues and proportions within and between individuals. Materials and MethodsUpper and lower first incisor to second premolar tooth rows were obtained from 30 individuals (n=300), from 3 archaeological samples. The teeth were micro-CT scanned and surface area and volumetric measurements were obtained from the surface meshes extracted. Dental variation of these measurements on a tooth and individual level was studied using Bayesian Multilevel Modelling. ResultsThe individual and tooth level variation differed by dental measurement, ranging between 9.5%-47.5% and 52.6-90.5% respectively. Enamel volume had the highest degree of individual-level variation in contrast to coronal dentine volume that had the lowest of individual-level variation. Tooth type, isomere, and position in field all showed a significant effect on the dental measurements examined in this study. DiscussionTooth selection and sampling strategies should consider individual and tooth-level variation, with at least one tooth from each type and isomere included in analyses. This will ensure that any population-level differences are not masked by variability between teeth. The low level of coronal dentine volume individual variation indicates that it is particularly useful in studies with small sample sizes.

Morphological covariation within the modern human postcanine dentition remains an open field of study. Analysis of covariation patterns of the three-dimensional (3D) shape between different tooth types has been seldom conducted, but it is relevant for the advancement of human biology and evolution, as well as dental anthropology, phylogeny, and medicine. Here, we analysed 3D shape covariation of the postcanine dentition (excluding third molars), both within and between dental arches using geometric morphometrics (GM). Based on high-resolution ({micro}CT) scans of 526 teeth from 136 individuals we found high pairwise correlation in tooth pairs within the dental arches (lower P3 and P4, r1 = 0.89; upper P3 and P4, r1 = 0.81; upper M1 and M2, r1 = 0.86). The correlation values between antagonists varied notably from the highest value detected between upper and lower M1s (r = 0.9), to the lowest between upper P4s and lower M1s (r = 0.58). Of all analysed tooth types, only the upper M1s showed moderate to high correlation in every pair analysis. Noticeably, unusually high covariation was detected between some of the tooth type pairs that do not articulate in a normal dentition (e.g., lower P3 and upper M2, r1 = 0.88). Furthermore, a relatively high covariation was found in the pairs of lower P4s and M1s (r1 = 0.79), and upper P4s and M1s (r1 = 0.77), which are the only tooth type pairs of the postcanine dentition belonging to different tooth classes (premolars and molars, respectively) and still serving similar masticatory functions. This study points to the fact that higher morphological integration seems to characterize teeth within the same dental arch rather than between antagonistic teeth. With this study, we provided an overview of pairwise correlations and strength of covariation between different tooth types. This information might inform future studies aimed at understanding developmental, phylogenetic, and functional aspects of the human postcanine dentition, including possible phenotype-genotype associations. However, with this study being the first one performed on a 3D sample of this size, we also report on obstacles and peculiarities that have been determined.

Focusing on craniofacial bones, this study investigates morphological variation related to sexual dimorphism in order to deepen our understanding of human biological diversity and to provide new data from a contemporary reference sample. Accordingly, the research was guided by three objectives (i) identify the facial regions exhibiting the greatest sexual dimorphism using landmark-based geometric morphometric method; (ii) evaluate the reliability of discriminant models based on these dimorphic regions; and (iii) conduct exploratory analyses to assign sex classification probabilities to ancient subjects using the discriminant models derived from a contemporary reference sample. The reference sample comprised 44 skulls from subjects of known sex who died in 2024. The ancient sample included 4 skulls recovered from the Grotte de La Medecine (France), attributed to the chalcolithic period. Fourteen facial landmarks were digitized using 3DSlicer. Generalized Procrustes Analysis was performed to extract shape variables and standardized coordinates for statistical analysis. Thin -Plate Spline transformations quantified and visualized deformation amplitudes between the female and male shapes. Landmarks in the orbital showed the highest deformation amplitudes. Goodall F-test comparing male and female shapes across three facial regions revealed significant sexual dimorphism only in the orbital region and the global facial shape. Discriminant analysis demonstrated that the orbital region provided the highest classification accuracy (88.6%) compared to the global facial region (72.7%). The discriminant models yielded high probabilities of male and female sex classification for the ancient subjects. Using a reproducible method, comparable levels of accuracy to those reported in the scientific literature were achieved despite a relatively modest sample size. These findings further confirm the significant role of the orbital region in human craniofacial sexual dimorphism.

Segmentation of high-resolution tomographic data is often an extremely time-consuming task and until recently, has usually relied upon researchers manually selecting materials of interest slice by slice. With the exponential rise in datasets being acquired, this is clearly not a sustainable workflow. In this paper, we apply the Trainable Weka Segmentation (a freely available plugin for the multiplatform program ImageJ) to typical datasets found in archaeological and evolutionary sciences. We demonstrate that Trainable Weka Segmentation can provide a fast and robust method for segmentation and is as effective as other leading-edge machine learning segmentation techniques.

In the past few decades, advances in 3D imaging have created new opportunities for reverse genetic screens. Rapidly growing datasets of 3D images of genetic knockouts require high-throughput, automated computational approaches for identifying and characterizing new phenotypes. However, exploratory, discovery-oriented image analysis pipelines used to discover these phenotypes can be difficult to validate because, by their nature, the expected outcome is not known a priori. Introducing known morphological variation through simulation can help distinguish between real phenotypic differences and random variation; elucidate the effects of sample size; and test the sensitivity and reproducibility of morphometric analyses. Here we present a novel approach for 3D morphological simulation that uses open-source, open-access tools available in 3D Slicer, SlicerMorph, and Advanced Normalization Tools in R (ANTsR). While we focus on diffusible-iodine contrast-enhanced micro-CT (diceCT) images, this approach can be used on any volumetric image. We then use our simulated datasets to test whether tensor-based morphometry (TBM) can recover our introduced differences; to test how effect size and sample size affect detectability; and to determine the reproducibility of our results. In our approach to morphological simulation, we first generate a simulated deformation based on a reference image and then propagate this deformation to subjects using inverse transforms obtained from the registration of subjects to the reference. This produces a new dataset with a shifted population mean while retaining individual variability because each sample deforms more or less based on how different or similar it is from the reference. TBM is a widely-used technique that statistically compares local volume differences associated with local deformations. Our results showed that TBM recovered our introduced morphological differences, but that detectability was dependent on the effect size, the sample size, and the region of interest (ROI) included in the analysis. Detectability of subtle phenotypes can be improved both by increasing the sample size and by limiting analyses to specific body regions. However, it is not always feasible to increase sample sizes in screens of essential genes. Therefore, methodical use of ROIs is a promising way to increase the power of TBM to detect subtle phenotypes. Generating known morphological variation through simulation has broad applicability in developmental, evolutionary, and biomedical morphometrics and is a useful way to distinguish between a failure to detect morphological difference and a true lack of morphological difference. Morphological simulation can also be applied to AI-based supervised learning to augment datasets and overcome dataset limitations.

A moving average smoothing method for extraction of cycles in time series data is described, with focus on obliquity cycles and fossil data. The proposed method is intended for cases where the environmental driver of phenotypic evolution can be shown to include obliquity cycles, either by power spectrum analysis or simply by inspection of raw or smoothed time series. The method gives improved mean trait predictions and better understanding when applied on stickleback fish fossil data from around 10 million years ago. The possibility to extract obliquity cycles will depend on the dynamics of the time series, and the method is thus not universally applicable. It may, however, be possible to adapt the size of the moving window to problems under study, or possibly to obtain improved predictions by inclusion of a sinusoidal component in the mean trait prediction modeling.

Understanding morphological variation is crucial for the study of speciation and for conservation as it helps in assessing biodiversity and predicting responses to environmental changes. These approaches are broadly applicable but are especially valuable in marine environments, where species are often elusive, difficult to study, and face heightened threats from rapid environmental shifts. The marine snail Littorina saxatilis is notable for its extensive polymorphism in shell shape, size, and colour, with ecotypes that evolve in response to environmental forces including wave exposure and crab predation. Morphometric tools have been central to investigating the mechanisms driving this phenotypic divergence; yet, a direct comparison of their methodological efficacy is lacking. Here, we took advantage of L. saxatilis ecotypes to contrast three morphometric approaches: elliptical Fourier analysis (EFA), landmarks-based geometric morphometrics (GM), and the growth-based model implemented in the ShellShaper software (SS). We assessed their clustering power, biological interpretability, robustness to measurement error and transferability among datasets. Our findings provide insights to guide method selection in studies aimed at exploring morphological variation: EFA is better suited for high-throughput screening and describing intermediate shapes; SS offers superior clustering power with highly interpretable growth parameters; and GM is best for detailed anatomical studies but is less efficient for large datasets. We provide guidelines to align method selection with specific research goals, balancing analytical efficiency with the required morphological and biological insight. By following this framework, researchers can ensure that robust morphological analysis is achieved, which is essential not only for elucidating mechanisms of adaptation and speciation but also for effective management and conservation of marine biodiversity.

AbstractAccounting for phylogenetic relatedness in the analysis of shape has become a common practice, deemed necessary to factor in the non-independence between species because of common ancestry. However, when adjusting error distributions to account for relatedness, the phylogenetic-generalised-least-squares (PGLS) test can obscure an important component of variation called conservative trait correlation (CTC). This is the amount of variation in a response variable that is both attributable to a predictor variable and phylogenetically structured. If CTC represents a large amount of correlated variation, true biological associations with strong phylogenetic signal (from unrepeated evolutionary events for example) might not be supported using a PGLS. We demonstrate this effect using geometric morphometric shape analysis on 370 crania from the speciose Australian rock- wallabies (genus Petrogale). In this clade, well-recognised allometric patterns such as scaling of the braincase (Hallers rule) and snout length (craniofacial evolutionary allometry) are supported using ordinary least squares (OLS) regression, but not PGLS, indicating that important between-species shape variation is lost. We then apply two methods capable of quantifying aspects of the missing variation: variation partitioning (VARPART), which estimates the proportion of variation shared between the predictor and phylogeny, and multi- response phylogenetic mixed models (MR-PMM), which identify the strength of correlation within the phylogenetic component of trait variance. Both methods show that CTC dominates the allometric shape variation in our sample, highlighting its importance in assessing phylogenetically informed models. We suggest approaches that can consider CTC become more widely used to better understand morphology and its predictors.

The general random walk model (GRW) of Hunt (2006) is used to infer directional evolution in mean trait values from sparse fossil data by modeling phenotypic change as the accumulated result of small steps with mean step sizes and step variances. Using simulations and real data cases, Ergon (2026) showed that the step variances can be estimated reasonably well only when the mean trait values have small measurement errors, while for fossil data with realistic measurement errors they appear to be extremely difficult to find, and they are often found to be negative. In the simulations Ergon (2026) assumed that the true phenotypic mean values were known. Here, I essentially repeat these simulations under the assumption that only mean trait values with large measurement errors are known, and based on weighted mean squared error (WMSE) comparisons the conclusion is that weighted least squares (WLS) is a better method than GRW. A second conclusion is that WLS is a better method also in the possibly rare cases with large measurement errors where the GRW parameters are estimated well. The GRW method is simply not flexible enough to handle such cases. A third conclusion is that Akaike Information Criterion (AIC) results for GRW models with large measurement errors relative to the step variance may be overly optimistic.

AO_SCPLOWBSTRACTC_SCPLOWThree dimensional morphometric methods are a powerful tool for comparative analysis of shape. However, morphological shape is often represented using landmarks selected by the user to describe features of perceived importance, and this may lead to over confident prediction of form-function relationships in subsequent analyses. We used Generalized Procrustes Analysis (GPA) of 13 homologous 3D landmarks and spherical harmonics (SPHARM) analysis, a homology-free method that describes the entire shape of a closed surface, to quantify the shape of the calcaneus, a landmark poor structure that is important in hind-limb mechanics, for 111 carnivoran species spanning 12 of 13 terrestrial families. Both approaches document qualitatively similar patterns of shape variation, including a dominant continuum from short/stout to long/narrow calcanea. However, while phylogenetic generalized linear models indicate that locomotor mode best explains shape from the GPA, the same analyses find that shape described by SPHARM is best predicted by foot posture and body mass without a role for locomotor mode, though effect sizes for all are small. User choices regarding morphometric methods can dramatically impact macroevolutionary interpretations of shape change in a single structure, an outcome that is likely exacerbated when readily landmarkable features are few.

Geometric morphometrics is widely employed across the biological sciences for the quantification of morphological traits. However, the scalability of these methods to large datasets is hampered by the requisite placement of landmarks, which can be laborious and time consuming if done manually. Additionally, the selected landmarks embody a particular hypothesis regarding the critical geometry pertinent to the biological inquiry at hand. Modifying this hypothesis lacks flexibility, necessitating the acquisition of an entirely new set of landmarks on the entire dataset to reflect any theoretical adjustments. In our research, we investigate the precision and accuracy of landmarks derived from the comprehensive set of functional correspondences acquired through the functional map framework of geometry processing. We use a deep functional map network to learn shape descriptors that effectively yield functional map-based and point-to-point correspondences between the specimens in our dataset. We then interrogate these maps to identify corresponding landmarks given manually placed landmarks from the entire dataset. We assess our method by automating the landmarking process on a dataset comprising mandibles from various rodent species, comparing its efficacy against MALPACA, a cutting-edge technique for automatic landmark placement. Compared to MALPACA, our model is notably faster and maintains competitive accuracy. The Root Mean Square Error (RMSE) analysis reveals that while MALPACA generally exhibits the lowest RMSE, our models perform comparably, especially with smaller training datasets, suggesting strong generalizability. Visual evaluations confirm the precision of our landmark placements, with deviations remaining within an acceptable range. These findings underscore the potential of unsupervised learning models in anatomical landmark placement, providing a viable and efficient alternative to traditional methods.

The Cooperative Eye Hypothesis (CEH) and the Gaze-Signalling Hypothesis (GSH) propose that the human eye--distinguished by pronounced scleral exposure and a uniformly white sclera-- evolved as a unique trait among primates to enhance eye-gaze visibility and facilitate cooperative communication. A recent review by Perea-Garcia and colleagues (2025) questioned four central premises of these hypotheses: (1) that human eye morphology is unique among primates, (2) that it is expressed consistently across individuals and populations, (3) that it improves gaze-following, and (4) that it is linked to the evolution of social cognition. Here, we revisit each claim through reevaluation of evidence and reanalysis of published data. First, we show that although some primates exhibit scleral depigmentation, humans uniquely combine this with high scleral exposure, resulting in markedly greater gaze visibility. Second, despite variation in scleral brightness, cross-species comparisons of sclera-iris-skin contrast confirm that human eyes remain distinct from those of other closely related apes. Third, experimental studies demonstrate that the white sclera confers a clear communicative advantage under ecologically relevant conditions, and that only humans consistently exploit these cues. Finally, developmental, cross-cultural, and neurocognitive evidence indicates that humans possess dedicated perceptual mechanisms for eye-gaze, consistent with its evolutionary embedding in social cognition. We conclude that while claims of human uniqueness should be moderated, the CEH and GSH remain the most plausible explanations for the evolution of human eye morphology. We also highlight key directions for future behavioral, anatomical, and genetic research.

Increasing variability down serially segmented structures, such as mammalian molar teeth and vertebrate limb segments, is a much-replicated pattern. The same phenotypic pattern has conflicting interpretations at different evolutionary scales. Macroevolutionary patterns are thought to reflect greater evolutionary potential in later-forming segments, but microevolutionary patterns are thought to reflect less evolutionary potential and greater phenotypic plasticity. We address this conflict by recalculating evolutionary potential (evolvability) from a systematic review of published mammalian molar sizes, then directly measure phenotypic plasticity from a controlled feeding experiment. Effects on lengths and widths are discordant in a way that suggests general growth pathways have a role in phenotypically plastic dental responses to nutrition. Effects on successive trait means do not necessarily increase downstream, contrary to long-standing hypotheses. We confirm prior findings of increasing non-inherited variance downstream, showing decoupling between effects on trait mean and variance. These patterns can be explained by a cascading model of tooth development compounding the effect of developmental instability as an influence separate from general environmental effects on the developing embryo. When evaluated in terms of evolvability, later-developing molars are equally or more evolvable than earlier-developing molars, aligning their microevolutionary potential with macroevolutionary patterns in other serially segmented structures.

Modern humans are acknowledged to have expanded across Earth from Africa. Biological measures have appeared to reflect this expansion, such as genetic diversity and cranial sexual size dimorphism. Admixture is known to be an issue for using diversity to locate where the expansion set out from in Africa; using cranial dimorphism should make admixture less of an issue. Therefore, cranial dimorphism could be of importance for clarifying the origin. This study used data sourced from the Howells dataset to infer if cranial form and shape dimorphisms indicate the expansion, and understand why cranial size dimorphism looks indicative. Form and shape dimorphisms were calculated through RMET. Size dimorphism was calculated beforehand. Cranial form and shape dimorphisms increased with distance from Africa. Form dimorphism suggested an area of origin which was predominantly in Africa and marginally in Asia. For shape dimorphism, locations spanned much farther beyond Africa. Hence, cranial form dimorphism seemed to be quite indicative of the expansion, unlike cranial shape dimorphism. Form is known to feature size and shape - cranial form dimorphism may signify the expansion mainly, or only, due to cranial size dimorphism. It seemed unsettled why cranial size and form dimorphisms seem related to the expansion because it was vague whether cranial size indicates the expansion more for males than females. Previously, a collective estimate of the origin (from biological measures including cranial size dimorphism) pointed to southern Africa; in the collective estimation process, cranial size dimorphism supported the south as would cranial form and shape dimorphisms.

The diversity of body shapes is one of the most prominent features of phenotypic variation in mammals. Yet, mammalian body shapes are poorly quantified and the underlying components contributing to its diversity as well as its relationship to other components of the skeleton are rarely tested. Here, we use lagomorphs (hares, rabbits and pikas) as a model system to (1) investigate which components of the skeleton contributed the most to body shape diversity, (2) examine the relationships between body shape and relative limb lengths, and (3) test how body size, ecotype, burrowing behavior, and locomotor mode influenced variation in lagomorph body shape and appendicular morphology. We quantified the body shape and functional proxies of the appendicular skeleton in 40 lagomorph species from osteological specimens held at museum collections. Using phylogenetic comparative methods, we found the relative length of the ribs and elongation or shortening of the thoracic and lumbar regions contributed the most to body shape evolution across lagomorphs. Second, we found that only leporids (hares and rabbits) exhibited a significant relationship between limb length and body shape, where more elongate species exhibit relatively shorter forelimbs and hindlimbs. Lastly, we found that models incorporating body size were the best predictors of lagomorph body shape and the majority of the appendicular traits, whereas models incorporating burrowing behavior and locomotor mode were largely poor fits. Broadly, these results indicate that larger lagomorphs tend to exhibit more robust body shapes with longer, more gracile forelimbs, whereas smaller lagomorphs tend to exhibit more elongate body shapes with shorter, more robust forelimbs. Overall, this work contributes to the growing understanding of mammalian body shape evolution and demonstrates the importance of not omitting body size in ecomorphological analyses.